In [55]:
import pandas as pd
import numpy as np
import os

import shap
import xgboost as xgb
from sklearn.linear_model import Ridge
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import precision_score, recall_score, accuracy_score, roc_curve, \
roc_auc_score, precision_recall_curve, confusion_matrix, r2_score, mean_absolute_error, \
mean_squared_error
import pandas_profiling as pp
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:60% !important; }</style>"))

Read the data

In [56]:
!dir data
 Volume in drive C is Windows
 Volume Serial Number is 5441-C8EF

 Directory of C:\Users\bialekj\Documents\Moje\Doktorat wdrozeniowy\wyklady\XAI\pd3\data

29.03.2020  18:24    <DIR>          .
29.03.2020  18:24    <DIR>          ..
22.03.2020  18:25         2˙452˙351 MEPS_data_preprocessed.csv
               1 File(s)      2˙452˙351 bytes
               2 Dir(s)  66˙544˙197˙632 bytes free
In [57]:
data_dir = 'data'
df = pd.read_csv(os.path.join(data_dir, 'MEPS_data_preprocessed.csv'))
df.reset_index(drop=True, inplace=True)
In [58]:
df_raw = df.copy()
In [59]:
df.head()
Out[59]:
PANEL REGION AGE31X GENDER RACE3 MARRY31X EDRECODE FTSTU31X ACTDTY31 HONRDC31 ... PCS42 MCS42 K6SUM42 PHQ242 EMPST31 POVCAT15 INSCOV15 INCOME_M HEALTHEXP PERSONWT
0 19 2 52 0.0 0.0 5 13 -1 2 2 ... 25.93 58.47 3 0 4 1 2 11390.0 46612 21854.981705
1 19 2 55 1.0 0.0 3 14 -1 2 2 ... 20.42 26.57 17 6 4 3 2 11390.0 9207 18169.604822
2 19 2 22 1.0 0.0 5 13 3 2 2 ... 53.12 50.33 7 0 1 2 2 18000.0 808 17191.832515
3 19 2 2 0.0 0.0 6 -1 -1 3 3 ... -1.00 -1.00 -1 -1 -1 2 2 385.0 2721 20261.485463
4 19 3 25 1.0 0.0 1 14 -1 2 2 ... 59.89 45.91 9 2 1 3 1 3700.0 1573 7620.222014

5 rows × 46 columns

Change categorical variables to strings (even though they're represented as numbers)

In [60]:
cats = ['REGION','MARRY31X','EDRECODE','FTSTU31X','ACTDTY31','HONRDC31',
            'RTHLTH31','MNHLTH31','HIBPDX','CHDDX','ANGIDX','MIDX','OHRTDX','STRKDX',
            'EMPHDX','CHBRON31','CHOLDX','CANCERDX','DIABDX','JTPAIN31','ARTHDX',
            'ARTHTYPE','ASTHDX','ADHDADDX','PREGNT31','WLKLIM31','ACTLIM31','SOCLIM31',
            'COGLIM31','DFHEAR42','DFSEE42','ADSMOK42','PHQ242','EMPST31','POVCAT15','INSCOV15']
In [61]:
df = df.copy()
for col in cats:
    df[col] = df[col].astype(str) 
In [62]:
# # Show summary report
# df.profile_report(title='Summary Report')

Let's have a closer look at the predicted variable (HEALTHEXP)

In [63]:
df['HEALTHEXP'].hist(bins=100)
Out[63]:
<matplotlib.axes._subplots.AxesSubplot at 0x16a52aaae48>
In [64]:
df['HEALTHEXP'].describe(percentiles=np.linspace(0,1,21))
Out[64]:
count     18350.000000
mean       5184.511608
std       15126.748532
min           0.000000
0%            0.000000
5%            0.000000
10%           0.000000
15%          20.350000
20%         114.000000
25%         198.000000
30%         302.000000
35%         431.000000
40%         585.000000
45%         779.000000
50%        1034.000000
55%        1339.000000
60%        1784.800000
65%        2343.850000
70%        3099.000000
75%        4219.500000
80%        5619.800000
85%        8137.000000
90%       12591.100000
95%       23417.000000
100%     659952.000000
max      659952.000000
Name: HEALTHEXP, dtype: float64
In [65]:
health_exp = df['HEALTHEXP'].values
In [66]:
plt.plot(np.linspace(0,100,21), np.percentile(health_exp, np.linspace(0,100,21)),'o');
plt.axhline(y=np.mean(health_exp), color='red', label = 'mean')
plt.xlabel('Percentile');
plt.ylabel('Healthexp');
plt.legend();
In [67]:
plt.plot(np.linspace(0,100,21), np.percentile(health_exp, np.linspace(0,100,21)),'o');
plt.axhline(y=np.mean(health_exp), color='red', label = 'mean')
plt.xlabel('Percentile');
plt.ylabel('Healthexp');
plt.yscale('log')
plt.legend();
In [68]:
plt.figure(figsize=(6,6))
df[['HEALTHEXP']].boxplot();
In [69]:
plt.figure(figsize=(6,6))
plt.boxplot(df['HEALTHEXP']);
plt.title('Health Exp')
plt.yscale('log')

Notes on predicted variable:

  • 14.3 % of observations are equal to 0.
  • The distribution is positevely skewed with many high-value outliers (very high treatment costs).
  • About 80 % of observations are below the mean. Std is three times higher than the mean.
In [70]:
cut_off = 100_000
In [71]:
df = df[df['HEALTHEXP']<cut_off]

Prepare features

In [72]:
# Drop panel number (not meant to be predictive) and sample weights
df.drop(columns = ['PANEL', 'PERSONWT'], inplace=True)
df.head()
Out[72]:
REGION AGE31X GENDER RACE3 MARRY31X EDRECODE FTSTU31X ACTDTY31 HONRDC31 RTHLTH31 ... ADSMOK42 PCS42 MCS42 K6SUM42 PHQ242 EMPST31 POVCAT15 INSCOV15 INCOME_M HEALTHEXP
0 2 52 0.0 0.0 5 13 -1 2 2 4 ... 2 25.93 58.47 3 0 4 1 2 11390.0 46612
1 2 55 1.0 0.0 3 14 -1 2 2 4 ... 2 20.42 26.57 17 6 4 3 2 11390.0 9207
2 2 22 1.0 0.0 5 13 3 2 2 1 ... 2 53.12 50.33 7 0 1 2 2 18000.0 808
3 2 2 0.0 0.0 6 -1 -1 3 3 1 ... -1 -1.00 -1.00 -1 -1 -1 2 2 385.0 2721
4 3 25 1.0 0.0 1 14 -1 2 2 1 ... 2 59.89 45.91 9 2 1 3 1 3700.0 1573

5 rows × 44 columns

In [73]:
y = df.pop('HEALTHEXP')

One-hot encoding for all the categorical features

In [74]:
cat_cols_one_hot = []
for col in cats:
    var_one_hot = pd.get_dummies(df[col],prefix=col, drop_first = True)
    df = df.drop(columns=[col])
    df = pd.concat([df, var_one_hot], axis=1)
    for cl in var_one_hot.columns:
        cat_cols_one_hot.append(cl)
In [75]:
df.head()
Out[75]:
AGE31X GENDER RACE3 PCS42 MCS42 K6SUM42 INCOME_M REGION_2 REGION_3 REGION_4 ... EMPST31_1 EMPST31_2 EMPST31_3 EMPST31_4 POVCAT15_2 POVCAT15_3 POVCAT15_4 POVCAT15_5 INSCOV15_2 INSCOV15_3
0 52 0.0 0.0 25.93 58.47 3 11390.0 1 0 0 ... 0 0 0 1 0 0 0 0 1 0
1 55 1.0 0.0 20.42 26.57 17 11390.0 1 0 0 ... 0 0 0 1 0 1 0 0 1 0
2 22 1.0 0.0 53.12 50.33 7 18000.0 1 0 0 ... 1 0 0 0 1 0 0 0 1 0
3 2 0.0 0.0 -1.00 -1.00 -1 385.0 1 0 0 ... 0 0 0 0 1 0 0 0 1 0
4 25 1.0 0.0 59.89 45.91 9 3700.0 0 1 0 ... 1 0 0 0 0 1 0 0 0 0

5 rows × 110 columns

More feature engineering in next HW ...

In [76]:
def train_print_save_models_reg(mods, dfTrain, dfTest, yTrain, yTest):
    trained_models = {}
    for model_name, mod_object in mods.items():
        arTrain = np.asarray(dfTrain)
        arTest = np.asarray(dfTest)
        md = mod_object.fit(arTrain,yTrain)

        y_pred_test = md.predict(arTest)
        y_pred_test[y_pred_test<0]=0
        y_pred_train = md.predict(arTrain)
        
        mae_train = mean_absolute_error(yTrain, y_pred_train)
        mae_test = mean_absolute_error(yTest, y_pred_test)
        rmse_train = np.sqrt(mean_squared_error(yTrain, y_pred_train))
        rmse_test = np.sqrt(mean_squared_error(yTest, y_pred_test))
        r2_train = r2_score(yTrain, y_pred_train)
        r2_test = r2_score(yTest, y_pred_test)
        print(model_name)
        print('Training R^2: {},  MAE:{}, RMSE:{} '.format(r2_train, mae_train, rmse_train))
        print('Test R^2:{}, MAE:{}, RMSE:{}'.format(r2_test, mae_test, rmse_test) )

        
        trained_models[model_name] = md
    return trained_models       
    
In [77]:
mods = {
    'XGB':xgb.XGBRegressor(max_depth=3),
    'LinReg':LinearRegression(),
    'Ridge_Reg':Ridge(alpha=.5)
}
In [78]:
dfTrain, dfTest, yTrain, yTest = train_test_split(df, y, random_state=42)
In [79]:
len(dfTrain), len(dfTest)
Out[79]:
(13713, 4571)
In [80]:
models = train_print_save_models_reg(mods, dfTrain, dfTest, yTrain, yTest)
XGB
Training R^2: 0.42262908558803425,  MAE:3920.1429059605084, RMSE:7798.7934522114665 
Test R^2:0.23218780618554702, MAE:4119.699536437045, RMSE:8135.991516875326
LinReg
Training R^2: 0.24761786167725786,  MAE:4445.026497647994, RMSE:8902.648071407939 
Test R^2:0.2522623347456887, MAE:4138.813566449553, RMSE:8028.928925100713
Ridge_Reg
Training R^2: 0.2476012849260607,  MAE:4445.127248945826, RMSE:8902.746144024033 
Test R^2:0.2524414329331225, MAE:4138.378630698163, RMSE:8027.967322656348
In [81]:
import tensorflow as tf
from keras.models import Sequential
from keras.optimizers import RMSprop
from keras.layers import Dense, Dropout, Embedding, Reshape
from keras.backend import clear_session
from keras.optimizers import Adam
from keras import backend as K
from keras.models import Model
from keras.layers import Dense
In [82]:
clear_session()
In [83]:
def r2(y_true, y_pred):
    SS_res =  K.sum(K.square( y_true-y_pred ))
    SS_tot = K.sum(K.square( y_true - K.mean(y_true) ) )
    return ( 1 - SS_res/(SS_tot + K.epsilon()) )
In [ ]:
 
In [87]:
model = Sequential([
    Dense(64, activation=tf.nn.relu, input_shape=[dfTrain.shape[1]]),
    Dropout(0.3),
    Dense(128, activation='relu'),
    Dense(64, activation='relu'),
    Dense(1)
  ])
In [88]:
optimizer = Adam(1e-4)
In [89]:
model.compile(loss='mae',
                optimizer = optimizer,
                metrics=[r2])
In [187]:
model.fit(x=dfTrain,
          y=yTrain, epochs=1000,
          batch_size =1024,
          validation_data=(dfTest,yTest),
          verbose=2,
          validation_freq=1, use_multiprocessing=True)

Simple ANN

In [91]:
y_pred = model.predict(dfTest)
In [92]:
r2_score(yTest, y_pred)
Out[92]:
0.12795497453172688
In [93]:
mae = mean_absolute_error(yTest, y_pred)
rmse = np.sqrt(mean_squared_error(yTest, y_pred))
print('MAE: {}, RMSE: {}'.format(mae, rmse))
MAE: 3526.694434853905, RMSE: 8670.665638462668
In [94]:
models['ANN'] = model

For some selected observation from this dataset, calculate the model predictions for model (1)

In [149]:
arTrain = np.asarray(dfTrain)
arTest = np.asarray(dfTest)
In [151]:
model1 = models['XGB']
model1
Out[151]:
XGBRegressor(base_score=0.5, booster=None, colsample_bylevel=1,
       colsample_bynode=1, colsample_bytree=1, gamma=0, gpu_id=-1,
       importance_type='gain', interaction_constraints=None,
       learning_rate=0.300000012, max_delta_step=0, max_depth=3,
       min_child_weight=1, missing=nan, monotone_constraints=None,
       n_estimators=100, n_jobs=0, num_parallel_tree=1,
       objective='reg:squarederror', random_state=0, reg_alpha=0,
       reg_lambda=1, scale_pos_weight=1, subsample=1, tree_method=None,
       validate_parameters=False, verbosity=None)
In [152]:
obs_num_1 = 10 
In [153]:
arobs = np.reshape(arTest[obs_num_1], (1,-1))
In [154]:
dfTest.iloc[obs_num_1,:]
Out[154]:
AGE31X        34.0
GENDER         0.0
RACE3          0.0
PCS42         -1.0
MCS42         -1.0
              ... 
POVCAT15_3     0.0
POVCAT15_4     0.0
POVCAT15_5     1.0
INSCOV15_2     0.0
INSCOV15_3     0.0
Name: 10429, Length: 110, dtype: float64
In [181]:
#wartość rzeczywista
yTest.iloc[obs_num_1]
Out[181]:
5833
In [183]:
#predykcja
model1.predict(arobs)[0]
Out[183]:
2412.7495

For an observation selected in (2), calculate the decomposition of model prediction using LIME / live / lime / localModel or similar technique (packages for R: live, live, localModel, iml, packages for python: lime).

In [156]:
import lime
from lime.lime_tabular import LimeTabularExplainer as LTE
In [157]:
cat_features = np.asarray([i for (i,col) in enumerate(dfTrain.columns) if col in cat_cols_one_hot])
In [158]:
feat_names = np.asarray(list(dfTrain.columns))
In [159]:
dfTrain.columns[cat_features]
Out[159]:
Index(['REGION_2', 'REGION_3', 'REGION_4', 'MARRY31X_10', 'MARRY31X_2',
       'MARRY31X_3', 'MARRY31X_4', 'MARRY31X_5', 'MARRY31X_6', 'MARRY31X_7',
       ...
       'EMPST31_1', 'EMPST31_2', 'EMPST31_3', 'EMPST31_4', 'POVCAT15_2',
       'POVCAT15_3', 'POVCAT15_4', 'POVCAT15_5', 'INSCOV15_2', 'INSCOV15_3'],
      dtype='object', length=103)
In [160]:
explainer = LTE(arTrain, mode='regression', feature_names=feat_names,
                class_names=['asd'], categorical_features=cat_features,
                verbose=True)
In [161]:
exp = explainer.explain_instance(arTest[10], model1.predict, num_features=6)
Intercept 31060.29384611636
Prediction_local [8464.22525601]
Right: 2412.7495
In [162]:
exp.show_in_notebook()

Opis zmiennych wybranych przez explainer i możliwe wartości:

  • PREGNT31 - (-1,1,2) czy w okresie referencyjnym jest/była w ciąży: -1 - nie dotyczy, 1 - tak, 2 - nie.
  • PCS42 - sumaryczna ocena zdrowia fizycznego (-1, wartość) - zmienna ciągła, jeśli wartość = -1 tzn., że ten wskaźnik nie mógł zostać obliczony: im wyżej tym lepiej.
  • DIABDX - cukrzyca (-1,1,2)
  • ACTLIM31 - ograniczenia związane z pracą/pracami domowymi (-1,1,2)
  • CANCERDX - nowotwór (-1,1,2)
  • RTHLTH31 - samoocena stanu zdrowia (-1, 1-5): im wyższa tym gorzej

Uwagi:

  • Model jako najistotniejszą cechę wybrał fakt bycia przez pacjenta w ciąży (w okresie, którego dotyczy zebrany wywiad). Wydaje się to mieć sens, jeśli kobieta była w ciąży to na pewno korzystała z usług medycznych - cały okres okołociążowy obiftuje w wizyty i badania lekarskie. Inna sprawa, że ta wybrana obserwacja dotyczy mężczyzny.

  • Druga wskazana cecha to fakt, że wskaźnik określający ogólny stan zdrowia fizycznego nie mógł zostać wyliczony (lub jest bardzo niski?). Być może zmienna ta wymaga modyfikacji.

  • Następnie model wskazuje na fakt niezdiagnozowania u pacjenta cukrzycy => wydatki na OM mniejsze.

  • Wyjaśnienie mówi też, że brak ograniczeń w wykonywaniu prac zarobkowych/domowych ujemnie wpływa na wydatki na opiekę medyczną.

  • Fakt niezdiagnozowania nowotworu ujemnie wpływa na przewidywane wydatki na OM.

  • Trudno jest zinterpretować wskazania modelu w przypadku samooceny stanu zdrowia (RTHLTH31). Wyjaśnienie mówi, że jeśli wartość ta jest różna od 4, to wydatki są generalnie mniejsze.

In [163]:
df_raw['RTHLTH31'].value_counts()
Out[163]:
 1    6060
 2    5699
 3    4290
 4    1688
 5     583
-1      30
Name: RTHLTH31, dtype: int64
In [164]:
dfTest.iloc[obs_num_1,:]
Out[164]:
AGE31X        34.0
GENDER         0.0
RACE3          0.0
PCS42         -1.0
MCS42         -1.0
              ... 
POVCAT15_3     0.0
POVCAT15_4     0.0
POVCAT15_5     1.0
INSCOV15_2     0.0
INSCOV15_3     0.0
Name: 10429, Length: 110, dtype: float64
In [165]:
df_raw['PHQ242'].hist(bins=2*7);
plt.plot()
Out[165]:
[]
In [166]:
several_obs = [23,32,42]
In [189]:
dfTest.iloc[several_obs,:]
Out[189]:
AGE31X GENDER RACE3 PCS42 MCS42 K6SUM42 INCOME_M REGION_2 REGION_3 REGION_4 ... EMPST31_1 EMPST31_2 EMPST31_3 EMPST31_4 POVCAT15_2 POVCAT15_3 POVCAT15_4 POVCAT15_5 INSCOV15_2 INSCOV15_3
9042 49 1.0 0.0 57.76 57.06 0 65000.0 0 1 0 ... 1 0 0 0 0 0 0 1 0 0
14462 81 1.0 0.0 39.34 59.42 0 8400.0 0 1 0 ... 0 0 0 1 0 0 0 0 0 0
12472 6 0.0 0.0 -1.00 -1.00 -1 0.0 0 1 0 ... 0 0 0 0 0 0 1 0 1 0

3 rows × 110 columns

In [167]:
for obs_id in several_obs:
    explainer.explain_instance(arTest[obs_id], model1.predict, num_features=6).show_in_notebook()
Intercept 27911.982779533977
Prediction_local [5443.54716479]
Right: 5787.1147
Intercept 29666.71042911351
Prediction_local [6688.25877657]
Right: 3926.0461
Intercept 31134.773394862666
Prediction_local [8382.62177426]
Right: 582.4159

Zmienne, które nie pojawiły się wcześniej:

  • MCS42 - sumaryczna ocena zdrowia psychicznego (-1, wartość) - zmienna ciągła, im wyżej tym lepiej.
  • WLKLIM31 - problemy z poruszaniem się (-1,1,2)
  • EDRECODE- wykształcenie (-1,1,2,13,14,15,16) - im wyżej tym lepsze.

Widzimy, że wśród tych kilku wybranych obserwacji pojawia się zestaw powtarzających się cech i ich wpływ na model (jakościowo) jest taki sam. Najistotniejsza pozostaje zmienna dotycząca ciąży. Wszystkie zmienne, które odpowiadają za fakt zdiagnozowania/niezdiagnozowania konkretnych chorób (kończące się na DX) wydają się sensownie wpływać na wyjścia modelu (jeśli ktoś ma/miał zdiagnozowaną daną chorobę to wydatki rosną). Dodatkowo pojawiają się zmienne związane z samopoczuciem i oceną stanu zdrowia, których interpretacja znów jest utrudniona.

Jeden z modeli wskazał, że fakt iż osoba na pewno nie jest doktorem/magistrem (EDRECODE_16=0) powoduje, że wydatki na OM są mniejsze.

Uwaga ogólna - kodowanie jedynkowe utrudnia odczyt takich wykresów.

In [168]:
model2 = models['ANN']
In [169]:
model2
Out[169]:
<keras.engine.sequential.Sequential at 0x16a510683c8>
In [170]:
ar = np.asarray(dfTest)
In [171]:
model2.predict(arobs)
Out[171]:
array([[258.94016]], dtype=float32)
In [172]:
exp1 = explainer.explain_instance(arTest[10], model1.predict, num_features=6)
exp2 = explainer.explain_instance(arTest[10], model2.predict, num_features=6)
Intercept 27109.911644629294
Prediction_local [7816.44345822]
Right: 2412.7495
Intercept 2451.989173279246
Prediction_local [1776.21085803]
Right: 258.94016
In [173]:
exp1.show_in_notebook()
exp2.show_in_notebook()

Zmienne, które nie pojawiły się wcześniej:

  • K6SUM52 - ogólna ocena odczuć (mentalnych) (-1, 0-24): im wyżej, tym gorzej.
  • Age31X - wiek
  • RACE3 - rasa - (0,1): 0 - white, 1 - black
In [186]:
df_raw['RACE3'].value_counts()
Out[186]:
0.0    12145
1.0     6205
Name: RACE3, dtype: int64
In [174]:
for obs_id in several_obs:
    explainer.explain_instance(arTest[obs_id], model2.predict, num_features=6).show_in_notebook()
Intercept 1579.647873352196
Prediction_local [1505.28022984]
Right: 1273.5063
Intercept 1314.1351448795226
Prediction_local [4060.30629061]
Right: 4641.412
Intercept 2483.839711174177
Prediction_local [1441.40882462]
Right: 276.10623

Porównanie modeli dla obserwacji nr 10 i kilku losowych:

  • Model ANN wskazał zupełnie inne cechy od modelu XGB
  • Model ANN wydaje się mniej stabilny - dla różnych obserwacji wskazuje zupełnie inne cechy. Widać więcej parametrów ciągłów, podczas gdy model XGB preferował parametry kategoryczne.